clear
clc
addpath(genpath(pwd));
%% User setting
dataOption        = 1;     %for using PMI 
startYear         = 1961;  %start year
hValue            = 1;     %forecast horizon
numBoot           = 5000;  %number of bootstraps
minObsRec         = 20;    %minimum recessions in bootstrapped sample
%% Loading data
loadData  = getData(dataOption,startYear);

% Recession probabilities
numResData     = sum(loadData.econ<0);
pHat           = sum(loadData.econ<0)/length(loadData.econ);
pHatSE         = sqrt(pHat*(1-pHat)/length(loadData.econ));
pHat_lower95   = pHat-1.654*pHatSE; %one-sided test
numRes_lower95 = pHat_lower95*length(loadData.econ);
pHat_lower99   = pHat-2.33*pHatSE;  %one-sided test
numRes_lower99 = pHat_lower99*length(loadData.econ);

% The maturies and yields
tmp       = loadData.yields_h1(1,:)<=120;
matSelect = loadData.yields_h1(1,tmp);
yields    = loadData.yields_h1(2:end,tmp);

% Selecting maturities
matSelectUsed = hValue:hValue:max(matSelect);
yieldsUsed    = zeros(size(yields,1),length(matSelectUsed));
for i=1:length(matSelectUsed)
    yieldsUsed(:,i) = yields(:,matSelectUsed(1,i) == matSelect);
end
numYields = size(yieldsUsed,2);

% We use a different short rate than implied by the GSW data set when
% hValue = 1 or hValue = 3
if hValue == 1
    % CRSP yield for h = 1
    yieldsUsed(:,1) = [loadData.crspRiskFree_h1(2:end,1);nan];
elseif hValue == 3
    % CRSP yield for h = 3
    yieldsUsed(:,1) = [loadData.crspRiskFree_h3(2:end,1);nan];
end

% Restrict data to a balanced panel
selectRow     = isnan(yieldsUsed)==0;
selectAll     = sum(selectRow,2)==size(selectRow,2);
yieldsUsed    = yieldsUsed(selectAll,:);
NBERrec       = loadData.NBERrec(selectAll);
econ          = loadData.econ(selectAll);
T             = size(yieldsUsed,1);
%% VAR model for bootstrap
[Wpca,Xpca] = pca(yieldsUsed,'NumComponents',3,'Centered',false);
resPCA      = yieldsUsed-Xpca*Wpca';
select      = isnan(resPCA);
sigmaError  = std(resPCA(select==0));
[B,OMEGA]   = Var1(Xpca);
% Bias-correcting the VAR
[h0Adj,hxAdj] = bootstrapVAR(Xpca',B(:,4),B(1:3,1:3),chol(OMEGA,'lower'),numBoot,1);
B = [hxAdj h0Adj];

% Equation for economic activity'
Yecon   = econ(2:end,1);
%Xecon   = [ones(T-1,1) Xpca(1:end-1,:) loadData.econ(1:end-1,1)];
Xecon   = [ones(T-1,1)  econ(1:end-1,1)];
results = nwest(Yecon,Xecon,hValue);

% The Var model: [PCA1;PCA2;PCA3;Z]
hx          = zeros(4,4);
hx(1:3,1:3) = B(1:3,1:3);
hx(4,4)     = results.beta(2);
h0          = [B(:,4); results.beta(1)];
nx          = 4;
factors     = [Xpca econ(1:end,1)]';
residuals   = zeros(size(factors));
for t=2:T
    residuals(:,t) =  factors(:,t) - (h0 + hx*factors(:,t-1));
end

%% The regression on empirical data
rec     = econ(1:end,1) < 0;
shortRateForXhr  = yieldsUsed(:,1);
shortRateOneYear = yieldsUsed(:,matSelectUsed==12);
resData = forecastRegressionYieldSpreadPCA_forBoot(yieldsUsed,shortRateForXhr,shortRateOneYear,matSelectUsed,rec,hValue);

% Autocorrelation of slope and slope*Rec
autoCorrSlope = autocorr(Xpca(:,2))
autoCorrSlopeRec = autocorr(Xpca(:,2).*rec)
autoCorrRec = autocorr(rec)
%% The bootstrap
XhrRegime_alfa_draws   = nan(numYields,numBoot);
XhrRegime_pca_draws    = nan(numYields,3,numBoot);
XhrRegime_dalfa_draws  = nan(numYields,numBoot);
XhrRegime_dbeta_draws  = nan(numYields,numBoot);
seXhrRegime_alfa_draws = nan(numYields,numBoot);
seXhrRegime_pca_draws  = nan(numYields,3,numBoot);
seXhrRegime_dalfa_draws= nan(numYields,numBoot);
seXhrRegime_dbeta_draws= nan(numYields,numBoot);
tXhrRegime_alfa_draws  = nan(numYields,numBoot);
tXhrRegime_pca_draws   = nan(numYields,3,numBoot);
tXhrRegime_dalfa_draws = nan(numYields,numBoot);
tXhrRegime_dbeta_draws = nan(numYields,numBoot);
b = 1;
idx = 0;
while b<=numBoot
   idx = idx +1;
   if mod(b,50) == 0 
       disp(['Draws so far: ',num2str(b)]);
   end
   % Bootstrapping the factors
   rng(idx,'twister')
   indexBoot          = randi(T-1,1,T)+1;
   factorsBoot        = zeros(nx,T);
   factorsBoot(:,1)   = factors(:,indexBoot(1,1));
   for s=2:T
       factorsBoot(:,s) = h0 + hx*factorsBoot(:,s-1) + residuals(:,indexBoot(1,s));
   end
   recBoot = real(factorsBoot(4,:)<0)';
   if sum(recBoot)>=minObsRec 
       % Yields in bootstrap
       yieldsBoot = factorsBoot(1:3,:)'*Wpca';
       yieldsBoot = yieldsBoot+repmat(sigmaError,size(yieldsBoot)).*randn(size(yieldsBoot));
       
       % The forecast regressions
       shortRateForXhrBoot  = yieldsBoot(:,1);
       shortRateOneYearBoot = yieldsBoot(:,matSelectUsed==12);
       resBoot = forecastRegressionYieldSpreadPCA_forBoot(yieldsBoot,shortRateForXhrBoot,shortRateOneYearBoot,matSelectUsed,recBoot,hValue);
              
       % Saving results
       XhrRegime_alfa_draws(:,b)   = resBoot.coef_alfa;
       XhrRegime_pca_draws(:,:,b)  = resBoot.coef_PCA';
       XhrRegime_dalfa_draws(:,b)  = resBoot.coef_dalfa;
       XhrRegime_dbeta_draws(:,b)  = resBoot.coef_dbeta;
       % The estimated SE
       seXhrRegime_alfa_draws(:,b)  = resBoot.se_alfa;
       seXhrRegime_pca_draws(:,:,b) = resBoot.se_PCA';
       seXhrRegime_dalfa_draws(:,b) = resBoot.se_dalfa;
       seXhrRegime_dbeta_draws(:,b) = resBoot.se_dbeta;
       % The t-stats
       tXhrRegime_alfa_draws(:,b)  = resBoot.tstat_alfa;
       tXhrRegime_pca_draws(:,:,b) = resBoot.tstat_PCA';
       tXhrRegime_dalfa_draws(:,b) = resBoot.tstat_dalfa;
       tXhrRegime_dbeta_draws(:,b) = resBoot.tstat_dbeta;
       
       b = b + 1;
   end
end

%% Analyzing the results
meanAlfaRegime   = mean(XhrRegime_alfa_draws,2); 
meanPCARegime    = squeeze(mean(XhrRegime_pca_draws,3));
seAlfaRegimeTrue = std(XhrRegime_alfa_draws,0,2); 
sePCARegimeTrue  = squeeze(std(XhrRegime_pca_draws,0,3));
tstatAlfaRegime  = prctile(tXhrRegime_alfa_draws,[2.5 97.5],2);
tstatPCARegime   = squeeze(prctile(tXhrRegime_pca_draws,[2.5 97.5],3));

meanAlfaRegime_change    = mean(XhrRegime_dalfa_draws,2); 
meanBetaRegime_change    = mean(XhrRegime_dbeta_draws,2);
seAlfaRegime_change_True = std(XhrRegime_dalfa_draws,0,2); 
seBetaRegime_change_True = std(XhrRegime_dbeta_draws,0,2);
tstatAlfaRegime_change   = prctile(tXhrRegime_dalfa_draws,[2.5 97.5],2);
tstatBetaRegime_change   = prctile(tXhrRegime_dbeta_draws, [2.5 97.5],2);

XhrRegime_alfa_draws_prctile  = prctile(XhrRegime_alfa_draws(:,:),[2.5 97.5],2);
XhrRegime_pca_draws_prctile   = prctile(squeeze(XhrRegime_pca_draws(:,:,:)),[2.5 97.5],3);
XhrRegime_dalfa_draws_prctile = prctile(XhrRegime_dalfa_draws(:,:),[2.5 97.5],2);
XhrRegime_dbeta_draws_prctile = prctile(XhrRegime_dbeta_draws(:,:),[2.5 97.5],2);

meanSEXhrRegime_alfa_draws    = mean(seXhrRegime_alfa_draws(:,:),2);
meanSEXhrRegime_pca_draws     = mean(seXhrRegime_pca_draws(:,:,:),3);
meanSEXhrRegime_dalfa_draws   = mean(seXhrRegime_dalfa_draws(:,:),2);
meanSEXhrRegime_dbeta_draws   = mean(seXhrRegime_dbeta_draws(:,:),2);
%% Ploting some results
close all
pcaPlot = 2;

figure('Name',['Point estimates: yield spread'],'NumberTitle','off')
subplot(2,2,1)
hold on
plot(matSelectUsed/12,meanAlfaRegime,'-k')
plot(matSelectUsed/12,XhrRegime_alfa_draws_prctile(:,1),':k')
plot(matSelectUsed/12,XhrRegime_alfa_draws_prctile(:,2),':k')
plot(matSelectUsed/12,resData.coef_alfa,'-r');
title('Regime-switching: Alfa')
axis tight
subplot(2,2,2)
hold on
plot(matSelectUsed/12,meanPCARegime(:,pcaPlot),'-k')
plot(matSelectUsed/12,XhrRegime_pca_draws_prctile(:,pcaPlot,1),':k')
plot(matSelectUsed/12,XhrRegime_pca_draws_prctile(:,pcaPlot,2),':k')
plot(matSelectUsed/12,resData.coef_PCA(pcaPlot,:),'-r');
title('Regime-switching: Beta')
axis tight
subplot(2,2,3)
hold on
plot(matSelectUsed/12,meanAlfaRegime_change,'-k')
plot(matSelectUsed/12,XhrRegime_dalfa_draws_prctile(:,1),':k')
plot(matSelectUsed/12,XhrRegime_dalfa_draws_prctile(:,2),':k')
plot(matSelectUsed/12,resData.coef_dalfa,'-r');
title('Regime-switching: Change in Alfa')
axis tight
subplot(2,2,4)
hold on
plot(matSelectUsed/12,meanBetaRegime_change,'-k')
plot(matSelectUsed/12,XhrRegime_dbeta_draws_prctile(:,1),':k')
plot(matSelectUsed/12,XhrRegime_dbeta_draws_prctile(:,2),':k')
plot(matSelectUsed/12,resData.coef_dbeta,'-r');
title('Regime-switching: Change in Beta')
axis tight


figure('Name',['SE: yields Spread'],'NumberTitle','off')
subplot(2,2,1)
hold on
plot(matSelectUsed/12,seAlfaRegimeTrue,'-k')
plot(matSelectUsed/12,meanSEXhrRegime_alfa_draws,'-r')
title('Regime-switching: Alfa')
legend('True SE','Mean of asymp. SE')
axis tight
subplot(2,2,2)
hold on
plot(matSelectUsed/12,sePCARegimeTrue(:,pcaPlot),'-k')
plot(matSelectUsed/12,meanSEXhrRegime_pca_draws(:,pcaPlot),'-r')
title(['Regime-switching: PCA',num2str(pcaPlot)])
axis tight
subplot(2,2,3)
hold on
plot(matSelectUsed/12,seAlfaRegime_change_True,'-k')
plot(matSelectUsed/12,meanSEXhrRegime_dalfa_draws,'-r')
title('Regime-switching: change in Alfa')
axis tight
subplot(2,2,4)
hold on
plot(matSelectUsed/12,seBetaRegime_change_True,'-k')
plot(matSelectUsed/12,meanSEXhrRegime_dbeta_draws,'-r')
title('Regime-switching: change in Beta')
axis tight


figure('Name',['t-stats: yield spread'],'NumberTitle','off')
subplot(2,2,1)
hold on
plot(matSelectUsed/12,tstatAlfaRegime(:,1),':k')
plot(matSelectUsed/12,tstatAlfaRegime(:,2),':k')
plot(matSelectUsed/12,resData.tstat_alfa,'-r')
title('Regime-switching: Alfa')
axis tight
subplot(2,2,2)
hold on
plot(matSelectUsed/12,tstatPCARegime(:,pcaPlot,1),':k')
plot(matSelectUsed/12,tstatPCARegime(:,pcaPlot,2),':k')
plot(matSelectUsed/12,resData.tstat_PCA(pcaPlot,:)','-r')
title(['Regime-switching: PCA',num2str(pcaPlot)])
axis tight
subplot(2,2,3)
hold on
plot(matSelectUsed/12,tstatAlfaRegime_change(:,1),':k')
plot(matSelectUsed/12,tstatAlfaRegime_change(:,2),':k')
plot(matSelectUsed/12,resData.tstat_dalfa,'-r')
title('Regime-switching: change in Alfa')
axis tight
subplot(2,2,4)
hold on
plot(matSelectUsed/12,tstatBetaRegime_change(:,1),':k')
plot(matSelectUsed/12,tstatBetaRegime_change(:,2),':k')
plot(matSelectUsed/12,resData.tstat_dbeta,'-r')
title('Regime-switching: change in Beta')

%% Saving results 
clear XhrRegime_alfa_draws XhrRegime_dalfa_draws XhrRegime_dbeta_draws XhrRegime_pca_draws
clear seXhrRegime_alfa_draws seXhrRegime_dalfa_draws seXhrRegime_dbeta_draws seXhrRegime_pca_draws
clear  tXhrRegime_alfa_draws  tXhrRegime_dalfa_draws tXhrRegime_dbeta_draws tXhrRegime_pca_draws

save(['GSWresPCA_numBoot',num2str(numBoot),'_minObsRec',num2str(minObsRec)])
